Import der Daten

library(tidyverse)
data_kidiq <- read_delim("https://raw.githubusercontent.com/sammerk/did_data/master/kidiq.csv", delim = ";")

Diejenigen, die nicht nutzen wollen, können den Datensatz hier im SPSS-Format herunterladen.

Parametrisierung

Zunächst soll der kid_score mit der metrischen Variable mom_iq prädiziert werden. In erfolgt das mit der Syntax

mod01 <- lm(kid_score ~ mom_iq, data = data_kidiq)
mod01

## 
## Call:
## lm(formula = kid_score ~ mom_iq, data = data_kidiq)
## 
## Coefficients:
## (Intercept)       mom_iq  
##       25.80         0.61

Grafisch kann dieses Modell wie folgt repräsentiert werden:

library(hrbrthemes)
data_kidiq %>% 
  ggplot(., aes(mom_iq, kid_score)) + 
  geom_point() + 
  stat_smooth(method = "lm", se = F) + 
  theme_ipsum()

## `geom_smooth()` using formula 'y ~ x'

Und die formale Beschreibung lautet:

\[ \operatorname{\widehat{kid\_score}} = 25.8 + 0.61(\operatorname{mom\_iq}) \]

Da der kid_score offensichtlich nicht einer Standardskalierung (z-Werte, IQ-Werte, t-Werte, …) entspricht, kann die Stärke des Effekts nur anhand einer Standardisierung bewertet werden.

mod02 <- lm(scale(kid_score) ~ scale(mom_iq), data = data_kidiq)
mod02

## 
## Call:
## lm(formula = scale(kid_score) ~ scale(mom_iq), data = data_kidiq)
## 
## Coefficients:
##   (Intercept)  scale(mom_iq)  
##    -5.068e-16      4.483e-01

Dies leisten auch packages wie {sjPlot} oder {stargazer} die darüber hinaus auch noch mehrere Modelle vergleichend darstellen können.

library(sjPlot)
tab_model(mod01, show.std = T, show.ci = F, show.p = F)

Nach Cohen (1988) gilt dabei:

• \(\beta \approx .1 \Rightarrow \text{kleiner Effekt}\)
• \(\beta \approx .3 \Rightarrow \text{moderater Effekt}\)
• \(\beta \approx .5 \Rightarrow \text{starker Effekt}\)

Inferenzstatistik

Unser \(\beta = .45\) beschreibt lediglich die Steigung einer Geraden, wenn man diese nach dem Kriterium der kleinsten Quadrate auf dem Datensatz kid_iq optimiert. Möchte man Schlussfolgerungen über den datengenerierenden Mechanismus anstellen benötigt man Inferenzstatistik.
p-Werte kann man etwa anhand der summary()-Funktion oder der tab_model() Funktion erhalten.

summary(mod01)

## 
## Call:
## lm(formula = kid_score ~ mom_iq, data = data_kidiq)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.753 -12.074   2.217  11.710  47.691 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.79978    5.91741    4.36 1.63e-05 ***
## mom_iq       0.60997    0.05852   10.42  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.27 on 432 degrees of freedom
## Multiple R-squared:  0.201,  Adjusted R-squared:  0.1991 
## F-statistic: 108.6 on 1 and 432 DF,  p-value: < 2.2e-16

Dabei sind zwei Inferenzstatistiken enthalten:

1. Jeder Koeffizient \(\beta_i\) wird gegen die 0 getestet.
2. Das Gesamtmodell wird gegen \(R^2 = 0\) getestet.

JZS-Bayes-Faktoren sind via {BayesFactor} ermittelbar:

library(BayesFactor)
lmBF(formula = kid_score ~ mom_iq, data = data_kidiq)

## Bayes factor analysis
## --------------
## [1] mom_iq : 4.516119e+19 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

Diagnostik der Voraussetzungen für die Inferenzstatistik

Diese Inferenzstatistiken sind nur unter Annahmen über den datengenerierenden Mechanismus berechenbar. Sind diese Annahmen verletzt sind die Inferenzstatistiken (mehr oder weniger stark) verzerrt. Daher nimmt die Diagnostik der Annahmen eine zentrale Bedeutung ein. Angenommen werden muss (mindestens)

• Linearität der Beziehung
• Normalverteilung der Residuen
• Homoskedastizität der Residuen
• Unabhängigkeit der Residuen

Eine gute Grundlage für die Diagnostik bietet das {performance}-package:

library(performance)
check_model(mod01)

Interpretation

Regressionsmodelle beschreiben immer nur bedingte Erwartungswerte von Datenpunkten (“für eine Gruppe von Müttern deren IQ im Schnitt X ist, ist ein kid_score von X am wahrscheinlichsten”). Woher diese Wahrscheinlichkeit rührt und wie gut mit dieser eine kausale Relationierung von X und Y gerechtfertigt werden kann hängt maßgeblich vom Design der Studie ab!

Parametrisierung

Die Regression ist insofern ein recht universelles Modellierungstool, als dass es auch die Aufnahme dichotomer Prädiktoren erlaubt. Bspw. kann man mit der Variable mom_hs der kid_score prädizieren:

ggplot(data_kidiq, aes(mom_hs, kid_score)) + 
  geom_point(alpha = .3) + 
  stat_smooth(method = "lm", se = F) + 
  theme_ipsum()

## `geom_smooth()` using formula 'y ~ x'

Die resultierenden Koeffizienten stellen die arithmetischen Mittelwerte der Gruppen dar und der p-Wert von \(\beta_1\) tested die \(H_0: EW(Gruppe_1) = EW(Gruppe_2)\).

mod03 <- lm(kid_score ~ mom_hs, data = data_kidiq)
mod03

## 
## Call:
## lm(formula = kid_score ~ mom_hs, data = data_kidiq)
## 
## Coefficients:
## (Intercept)       mom_hs  
##       77.55        11.77

Wobei \[ \operatorname{kid\_score} = \beta_{0} + \beta_{1}(\operatorname{mom\_hs}) + \epsilon \]

Diagnostik der Voraussetzungen für die Inferenzstatistik

Da Regressionsmodelle mit Dummyvariablen eben auch Regressionsmodelle sind 😃 gilt es auch bei diesen dieselben Voraussetzungen zu prüfen:

check_model(mod03)

Parametrisierung

In Regressionsmodelle können problemlos mehrere Prädiktoren aufgenommen werden.

mod04 <- lm(kid_score ~ mom_iq + mom_hs, data = data_kidiq)
summary(mod04)

## 
## Call:
## lm(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.873 -12.663   2.404  11.356  49.545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.73154    5.87521   4.380 1.49e-05 ***
## mom_iq       0.56391    0.06057   9.309  < 2e-16 ***
## mom_hs       5.95012    2.21181   2.690  0.00742 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.14 on 431 degrees of freedom
## Multiple R-squared:  0.2141, Adjusted R-squared:  0.2105 
## F-statistic: 58.72 on 2 and 431 DF,  p-value: < 2.2e-16

Grafisch kann dies dann im entsprechenden n-dimensionalen Raum repräsentiert werden (siehe Folien aus Universaltool Regressionsanalyse I). In der formalen Repräsentation wird schlicht das Polynom um einen Summanden erweitert: \[ \operatorname{kid\_score} = \beta_{0} + \beta_{1}(\operatorname{mom\_iq}) + \beta_{2}(\operatorname{mom\_hs}) + \epsilon \]

Interpretation

Besonders interessant an multiplen Regressionsmodellen ist, dass bestimmte kausale Mechanismen falsifiziert werden können. Dazu werden dann meist die Ergebnisse der einfachen Modelle (mit jeweils einem Prädiktor) mit den Parameterschätzungen im multiplen Modell verglichen:

tab_model(mod01, mod03, mod04, show.ci = F)

So können Annahmen über Supressor- oder Konfounderkonstellationen falsifiziert werden. Hier etwa wird klar, dass die Daten mit dem Mechanismus A nicht im Einklang stehen, aber mit dem Mechanismus B und C.

Diagnostik der Voraussetzungen für die Inferenzstatistik

Auch Modelle der multiplen Regression unterliegen den gleichen 4 Annahmen. Hinzu kommt (je nach Methode der Schätzung), dass keine Kolinearität vorliegen sollte. check_model() nimmt dies Annahme direkt in die Prüfung mit auf, wenn das Argument der Funktion ein Modell mit mehreren Prädiktoren darstellt.

check_model(mod04)

Typischerweise werden dann wieder p-Werte für die Hypothesen \(\beta_i = 0\) und $ R^2 = 0$ berechnet.
Bayes Faktoren beziehen sich jedoch immer auf das inkrementelle \(R^2\), also die Frage: Klärt ein Prädiktor oder eine Gruppe von Prädiktoren zuätzliche Varianz auf?

# Test gegen intercept only modell
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)

## Warning: data coerced from tibble to data frame

## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 1.508721e+20 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

# Test gegen Modell mit mom_hs als Prädiktor
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)/
  lmBF(formula = kid_score ~ mom_hs, data = data_kidiq)

## Warning: data coerced from tibble to data frame

## Warning: data coerced from tibble to data frame

## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 7.578673e+15 ±0%
## 
## Against denominator:
##   kid_score ~ mom_hs 
## ---
## Bayes factor type: BFlinearModel, JZS

# Test gegen Modell mit mom_iq als Prädiktor
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)/
  lmBF(formula = kid_score ~ mom_iq, data = data_kidiq)

## Warning: data coerced from tibble to data frame

## Warning: data coerced from tibble to data frame

## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 3.340747 ±0%
## 
## Against denominator:
##   kid_score ~ mom_iq 
## ---
## Bayes factor type: BFlinearModel, JZS

Parametrisierung

Interessiert man sich dafür, wie sich die Effekte je nach Ausprägung einer weiteren nominalen Variablen unterscheiden, kann man Modelle mit Interaktionen eines metrischen und eines dichotomen Prädiktors spezifizieren.

ggplot(data_kidiq, aes(mom_iq, kid_score, 
                       color = as.factor(mom_hs), 
                       group = as.factor(mom_hs))) + 
  geom_point() + 
  stat_smooth(method = "lm") + 
  theme_ipsum()

## `geom_smooth()` using formula 'y ~ x'

Die entsprechenden Koeffizienten des Modells mit Interaktionseffekt sind jedoch wesentlich schwieriger zu interpretieren:

mod05 <- lm(kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = data_kidiq)
mod05

## 
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = data_kidiq)
## 
## Coefficients:
##   (Intercept)         mom_hs         mom_iq  mom_hs:mom_iq  
##      -11.4820        51.2682         0.9689        -0.4843

\[ \operatorname{kid\_score} = \alpha + \beta_{1}(\operatorname{mom\_hs}) + \beta_{2}(\operatorname{mom\_iq}) + \beta_{3}(\operatorname{mom\_hs} \times \operatorname{mom\_iq}) + \epsilon \]

Leichter wird dies, wenn man den kontinuierlichen Prädiktor zentriert oder standardisiert.

library(equatiomatic)
data_kidiq <- data_kidiq %>% 
  mutate(mom_iq_centered = as.numeric(scale(mom_iq, center = T, scale = F)),
         mom_iq_zstand = as.numeric(scale(mom_iq, center = T, scale = T)),
         mom_age_zstand = as.numeric(scale(mom_age, center = T, scale = T)))

mod06 <- lm(kid_score ~ mom_hs + mom_iq_zstand + mom_hs:mom_iq_zstand, data = data_kidiq)

mod06

## 
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq_zstand + mom_hs:mom_iq_zstand, 
##     data = data_kidiq)
## 
## Coefficients:
##          (Intercept)                mom_hs         mom_iq_zstand  
##               85.407                 2.841                14.533  
## mom_hs:mom_iq_zstand  
##               -7.264

ggplot(data_kidiq, aes(mom_iq_zstand, kid_score, 
                       color = as.factor(mom_hs), 
                       group = as.factor(mom_hs))) + 
  geom_point() + 
  stat_smooth(method = "lm") + 
  theme_ipsum()

## `geom_smooth()` using formula 'y ~ x'

extract_eq(mod06, use_coefs = T, terms_per_line = 2, wrap = T)

\[ \begin{aligned} \operatorname{\widehat{kid\_score}} &= 85.41 + 2.84(\operatorname{mom\_hs})\ + \\ &\quad 14.53(\operatorname{mom\_iq\_zstand}) - 7.26(\operatorname{mom\_hs} \times \operatorname{mom\_iq\_zstand}) \end{aligned} \]

check_model(mod06)